*This do file creates figured from stored estimate files and aggregated datasets that were created on IRS server

*net from "https://mloeffler.github.io/stata/"
*net install wordwrap

clear all
set more off
set matsize 11000
set emptycells drop
set scheme s2mono	


*The directory path has been replaced with "****"
global path "****/Results"
global output "****/Figures"

cd "$path"

/* =====================
Figure I: Coverage Effect by Age 
========================*/

clear all
mat drop _all
estimates clear

forval c = 2/7{
	estimates use "$path/Estimates/cov1718_age10_033020.ster", number(`c')
	mat estimates = nullmat(estimates) \ (_b[treatment], (_b[treatment] - 1.96*_se[treatment]), (_b[treatment] + 1.96*_se[treatment]))
}

set obs 6 
gen age_group = _n
svmat estimates

rename estimates1 effect 
rename estimates2 lower
rename estimates3 upper

local note = "{it:Notes}: The figure displays the estimated effect of the intervention on coverage by age. The outcome is the number of months of coverage enrolled in during 2017-18. Each specification is limited to individuals whose ages fell into the specified range at the end of 2017. The figure excludes individuals with full coverage in January through November of 2016. Brackets denote the 95% confidence interval based on standard errors that are clustered by household."
wordwrap "`note'", l(90)	

twoway 	   (scatter effect age_group , clpattern(1) clcolor(gs1) msymbol(o) mcolor(black)) ///
		   (rcap upper lower age_group, clpattern(1) lcolor(black) lwidth(thin)), ///
			xlabel(1 `" 5-14  "'  4 `" 35-44 "' ///
				   2 `" 15-24 "'  5 `" 45-54 "' ///
				   3 `" 25-34 "'  6 `" 55-64 "', labsize(small)) ///			
			xtitle("Age Group", height(4)) ///
			ytitle("Effect of Intervention on Coverage (months)", size(medsmall) height(5)) ///
			graphregion(fcolor(white) color(white) lcolor(white)) ///
	        legend(off) ///
			note("`r(text)'")

graph export "$output/Coverage Effect by Age (Intensive Margin).eps", replace as(eps) 


/* ===============================
Figure II: Effect of Intervention on Coverage
==================================*/

* Panel A

use "$path/Data/monthly_cov_notall.dta", clear 

* Aggregate by control/treat for each type of coverage

forval x = 1/36 {
	rename month`x'_se month_se`x'
}	

* Reshape long
reshape long month month_se, i(treatment) j(num)
rename month cov
rename month_se cov_se
rename num month

* Change to percentage points
replace cov = cov*100

* Figure
twoway (connected cov month if treatment == 0, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)) ///
	   (connected cov month if treatment == 1, msymbol(d) mcolor(gs1) clcolor(gs1)) , ///
		ysc(range(10 55)) ylabel(10(5)50, format(%5.0f) labsize(small)) /// 
		xline(13, lpattern(shortdash) lcolor(black) lwidth(medium)) ///
		xlabel(1 `""Jan" "2016" "'  					 ///
				5"May"  9"Sep"   	 ///
				13 `""Jan" "2017" "'					 ///
				17"May" 21"Sep"  ///
				25 `""Jan" "2018" "'					 ///
				29"May"  33"Sep"   ///
				,  format(%tmMon) labsize(small))  ///
		xtitle("", height(4)) ///
		ytitle("Coverage Rate (%)", size(medsmall)) ///
		title("Panel A", size(medsmall) height(0)) ///
		legend(label(1 "Control") label(2 "Treatment") rows(2) pos(0) bplace(neast) order(2 1)) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		saving(1a, replace)

* Panel B
		
mat drop _all
forval c = 1/36 {
estimates use "$path/Estimates/monthly_cov_effect_051720.ster", number(`c')
mat results = nullmat(results) \ (_b[treatment]*100, 100*(_b[treatment] - 1.96*_se[treatment]) ,100*(_b[treatment] + 1.96*_se[treatment]))
}

clear
set obs 36
gen month = _n
svmat results

rename results1 diff
rename results2 lower 
rename results3 upper

* Figure
twoway (connected diff month, msymbol(o) mcolor(gs1) clcolor(gs1)) ///
	   (rcap upper lower month,  lcolor(gs1) lwidth(vthin)), ///
	     ysc(range(0 2.2)) ylabel(0(0.4)2.0, format(%5.1f) labsize(small)) ///
		xline(13, lpattern(shortdash) lcolor(black) lwidth(medium)) ///
		xlabel(1 `""Jan" "2016" "'  					 ///
				5"May"  9"Sep"   	 ///
				13 `""Jan" "2017" "'					 ///
				17"May" 21"Sep"  ///
				25 `""Jan" "2018" "'					 ///
				29"May"  33"Sep"   ///
				,  format(%tmMon) labsize(small))  ///
		xtitle("", height(4)) ///
		ytitle("Difference in Coverage Rate (p.p.)", size(medsmall)) ///
		title("Panel B", size(medsmall) height(0)) ///
		legend(label(1 "Difference") label(2 "95% CI") ///
					   subtitle("") rows(2) pos(0) bplace(neast)) ///	
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		saving(1b, replace)

		
* Combine Panel A and B
* ----------------------
graph combine 1a.gph 1b.gph, col(1) ysize(9) xsize(7) graphregion(margin(l=5 r=5)) ///
	note("{it:Notes}: Panel A displays the shares of the treatment and control group enrolled in any coverage" ///
			 "in the specified month. Panel B displays the difference in the share between the treatment and" ///
			 "control groups enrolled in any coverage in the specified month. Units are percentage points" ///
			 "(0-100). Both panels exclude individuals with full coverage in January through November of" ///
			 "2016. Brackets denote the 95% confidence interval based on standard errors that are clustered" ///
			 "by household.", ///
					size(vsmall)  height(15) linegap(0.5)) ///
					graphregion(fcolor(white) color(white) lcolor(white))
					
graph export "$output/Coverage Effect Over Time.eps", replace as(eps) 
cap erase 1a.gph 
cap erase 1b.gph 


			
/* ===========================================
Figure III: Effect of Intervention on Middle Age Mortality
==============================================*/

* Panel A:

use "$path/Data/mortality_by_month_051920.dta", clear
	
	*Sample Restriction
	keep if notall16 == 1

	* Calculate rates
	forval x = 1/36 {
		replace dead_by_`x' = (dead_by_`x' / dep_tin)*100
	}
	
	* Reshape
	reshape long dead_by_, i(treatment) j(month) 
	rename dead_by_ percent_dead
	format dep_tin %12.0f
	rename dep_tin count
	
	* Figure
	twoway (connected percent_dead month if treatment == 1, clpattern(1) clcolor(gs1) msymbol(o) msize(small) mcolor(gs1)) ///
		   (connected percent_dead month if treatment == 0, clpattern(dash) clcolor(black)  msymbol(oh) msize(small) mcolor(gs1)), ///
			xline(13, lpattern(shortdash) lcolor(black) lwidth(medium)) ///
			ysc(range(0 1.5)) ylabel(0(0.2)1.4, format(%5.1f) labsize(small)) ///
			xlabel(1 `" "Jan" "2016" "'  3 `" "Mar" "' ///
				   5 `" "May"  "' 		 7 `" "Jul"  "' ///
				   9 `" "Sep"  "' 		 11 `" "Nov"  "' ///
				   13 `" "Jan" "2017" "' 15 `" "Mar"  "' ///
				   17 `" "May"  "' 		 19 `" "Jul"  "' ///
				   21 `" "Sep" "' 		 23 `" "Nov"  "' ///
				   25 `" "Jan" "2018" "' 27 `" "Mar"  "' ///
				   29 `" "May"  "'		 31 `" "Jul"  "' ///
				   33 `" "Sep" "' 		 35 `" "Nov"  "', ///
					labsize(small)) ///			
			xtitle("", height(0)) ///	
			title("Panel A", size(medsmall) height(0)) ///
			ytitle("Cumulative Mortality (%)", size(medsmall) height(4)) ///
			legend(label(1 "Treatment") label(2 "Control") ///
						   subtitle("") size(small) rows(2) pos(0) bplace(nwest)) ///	
			graphregion(fcolor(white)) /// 
			saving(a, replace) 

		
			
* Panel B:

estimates clear
eststo clear

cap mat drop estimates
forval c = 1/7{
estimates use "$path/Estimates/mort_by_month_051920.ster", number(`c')
mat estimates = nullmat(estimates) \ (_b[treatment]*100, 100*(_b[treatment] - 1.96*_se[treatment]) ,100*(_b[treatment] + 1.96*_se[treatment]))
}

clear
set obs 7
gen month = (_n-1)*6 + 1 if _n < 7
replace month = 36 if _n == 7
svmat estimates

rename estimates1  percent_dead
rename estimates2 upper
rename estimates3 lower


	* Figure
	twoway (connected percent_dead month , clpattern(1) clcolor(gs1) msymbol(o) mcolor(black)) ///
		   (rcap upper lower month, clpattern(1) lcolor(black) lwidth(thin)), ///
			xline(13, lpattern(shortdash) lcolor(black) lwidth(medium)) ///
			yline(0, lwidth(thin) lcolor(black)) ///
			ysc(range(-0.15 0.07)) ylabel(-0.15(0.05)0.05, format(%5.2f) labsize(small)) ///
			xlabel(1  `" "Jan" "2016" "'  7  `" "Jul" "' ///
				   13 `" "Jan" "2017" "'  19 `" "Jul"  "' ///
				   25 `" "Jan" "2018" "'  31 `" "Jul"  "' ///
				   36 `" "Dec" "' , ///
					labsize(small)) ///	
			xtitle("", height(0)) ///
			ytitle("Difference in Cumulative Mortality (p.p.)", size(medsmall) height(4)) ///
			title("Panel B", size(medsmall) height(0)) /// 
			legend(label(1 "Difference") label(2 "95% CI") ///
						   subtitle("") size(small) rows(2) pos(0) bplace(nwest) )  ///
			graphregion(fcolor(white)) ///
			saving(b, replace)

local note = "{it:Notes}: Panel A displays the share of middle age adults that died during or prior to the specified month. Panel B displays the difference in the cumulative mortality rate among middle age adults in the control and treatment groups. The difference is calculated at six-month intervals that extend through the end of the specified month. Units are percentage points (0-100). Both panels are limited to individuals between the ages of 45-64 at the end of 2017 and exclude individuals with full coverage in January through November of 2016. Brackets denote the 95% confidence interval based on standard errors that are clustered by household." 
wordwrap "`note'", l(90)			

* Combine Panel A and B
* ----------------------
graph combine a.gph b.gph, col(1) ysize(10) xsize(8) graphregion(margin(l=5 r=5))  ///
	note("`r(text)'", ///
		  linegap(0.5) size(vsmall)) ///
		 graphregion(fcolor(white) color(white) lcolor(white)) 	 
	
graph export "$output/Mortality Over Time.eps", replace as(eps)

cap erase a.gph 
cap erase b.gph 


/* ============================================
Figure IV: Coverage and Mortality Effects by Treatment Arm and Individual Characteristics
=============================================== */

use "$path/Data/arms_collapse_050520.dta", clear

rename dep_tin count
drop any_covered2017 covered2017 any_covered1718
egen treatment_arm = group(base_no_sp early_no_sp nonpersonalize_no_sp exemption_info_no_sp base_sp early_sp nonpersonalize_sp exemption_info_sp)
collapse (rawsum) count (mean) covered1718 dead1718 [fw = count], by(age_45_64 above_med_cov_effect cov2016_11 treatment_arm)

preserve
keep if treatment_arm == 1
rename count count_c
rename covered1718 covered1718_c
rename dead1718 dead1718_c
tempfile control_group
save "`control_group'"
restore 

keep if treatment_arm != 1
merge m:1 above_med_cov_effect cov2016_11 age_45_64 using "`control_group'", nogen

foreach i of var covered1718 dead1718{
gen effect_`i' = `i' - `i'_c 
}

replace effect_dead1718 = effect_dead1718*100

egen sum = sum(count)

gen weight = count/sum 


reg effect_dead1718 effect_covered1718 [pw = weight], robust


local note = "{it:Notes}: The figure plots the estimated effect of the intervention on coverage and mortality, separately for groups of individuals defined based on treatment arm, state of residence, age, and prior-year coverage. With respect to state of residence, individuals are grouped according to whether the estimated coverage effect for their state was above or below the median state's estimated coverage effect. With respect to age, individuals are grouped according to whether they are aged 45-64 at the end of 2017; individuals aged older than 64 are excluded from the analysis. With respect to prior year coverage, individuals are grouped according to whether they were enrolled in full coverage in January through November of 2016. The coverage effect (x-axis) corresponds to the effect of the intervention on the number of months of coverage enrolled in during 2017-18. The mortality effect (y-axis) corresponds to the effect of the intervention on mortality in 2017-18; units are percentage points (0-100). The best linear fit (dashed line), weighted by cell sample size, has a slope of  -0.081, with standard error 0.022."
wordwrap "`note'", l(90)

twoway 	 (lfit effect_dead1718 effect_covered1718 [pw = weight], lpattern(dash) lcolor(black) lwidth(thin)) ///
		 (scatter effect_dead1718 effect_covered1718 [fw = count], ///
			clpattern(1) clcolor(gs1) msize(tiny) msymbol(circle) mcolor(black)), ///
			graphregion(fcolor(white) color(white) lcolor(white)) ///
			xtitle("Coverage Effect (Months)", height(4)) ///
			xsc(r(0 0.6)) xlabel(0(0.2)0.6) ///
			ytitle("Mortality Effect (p.p.)", height(5)) ///
			legend(off) ///
			note("`r(text)'")

graph export "$output/Effect by Covariate and Treatment Arms.eps", replace as(eps) 



/* ==================================================
Figure V: Effect of Coverage on Mortality by Age
===================================================*/

clear all
mat drop _all
estimates clear

forval c = 2/7{
	estimates use "$path/Estimates/iv1718_age10_033020.ster", number(`c')
	mat estimates = nullmat(estimates) \ (_b[covered1718]*100, 100*(_b[covered1718] - 1.96*_se[covered1718]) ,100*(_b[covered1718] + 1.96*_se[covered1718]))
}

set obs 6 
gen age_group = _n
svmat estimates

rename estimates1 effect 
rename estimates2 lower
rename estimates3 upper

local note = "{it:Notes}: The figure displays the estimated effect of coverage on mortality by age. Each estimate is obtained from an IV specification in which an indicator for treatment group assignment is used to instrument for months of coverage during 2017-18. The outcome is mortality in 2017-18; units are percentage points (0-100). Each specification is limited to individuals whose ages fell into the specified range at the end of 2017. The figure excludes individuals with full coverage in January through November of 2016. Brackets denote the 95% confidence interval based on standard errors that are clustered by household." 
wordwrap "`note'", l(88)	

twoway 	   (scatter effect age_group , clpattern(1) clcolor(gs1) msymbol(o) mcolor(black)) ///
		   (rcap upper lower age_group, clpattern(1) lcolor(black) lwidth(thin)), ///
			yline(0, lwidth(thin) lcolor(black)) ///		
			xlabel(1 `" 5-14  "'  4 `" 35-44 "' ///
				   2 `" 15-24 "'  5 `" 45-54 "' ///
				   3 `" 25-34 "'  6 `" 55-64 "', labsize(small)) ///			
			xtitle("Age Group", height(4)) ///
			ytitle("Effect of Coverage on Mortality (p.p.)", size(medsmall) height(5)) ///
			graphregion(fcolor(white) color(white) lcolor(white)) ///
	        legend(off) ///
			note("`r(text)'")


graph export "$output/2SLS Effect by Age.eps", replace as(eps) 
						

*****************************Appendix Figures***********************************


/*===================
Figure A.I: Social Security Death Index Data Check
=====================*/
clear all
import excel "$path/Data/Mortality vital stats vs DM1 (raw).xlsx", sheet("Sheet1") firstrow 

rename C year_month
drop if year_month < 201001
gen obs_no = _n

replace VitalStats = VitalStats/100000
replace DM1 = DM1/100000


local note = "{it:Notes}: The figure compares monthly deaths recorded in the Social Security Death Index and the National Vital Statistics System Mortality Data from 2010 through the end of 2018. Deaths that occur abroad are recorded in the Social Security Death Index but not in the National Vital Statistics System."
wordwrap "`note'", l(90)	

twoway (line VitalStats obs_no, lpattern(shortdash) lcolor(gs1) lwidth(vthin)) ///
	   (line DM1 obs_no, lpattern(solid) lcolor(gs1) lwidth(vthin)),  ///	
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		xlabel(1 "2010" 13 "2011"  25 "2012" 37 "2013"	///
			   49 "2014" 61 "2015" 73 "2016" 85 "2017"	///
			   97 "2018" ,  format(%tmMon) labsize(small)) ///
		ytitle("Deaths (in 100,000's)", size(medsmall) height(6)) ///
		xtitle("", height(4)) ///
		ysc(range(0 3)) ylabel(0(0.5)3) ///
		legend(label(1 "National Vital" "Statistics System") label(2 "Social Security Death Index") ///
				rows(2) size(vsmall) pos(0) bplace(seast))  ///
		note("`r(text)'")

graph export "$output/Social Security Death Index Check.eps", replace as(eps) 


			

/* ===============================
Figure A.III: Coverage Effect Among Middle Age Adults
==================================*/


*Panel A

use "$path/Data/monthly_cov_notall_4564.dta", clear 

forval x = 1/36 {
	rename month`x'_se month_se`x'
}	

* Reshape long
reshape long month month_se, i(treatment) j(num)
rename month cov
rename month_se cov_se
rename num month

* Change to percentage points
replace cov = cov*100

* Figure
twoway (connected cov month if treatment == 0, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)) ///
	   (connected cov month if treatment == 1, msymbol(d) mcolor(gs1) clcolor(gs1)) , ///
		ysc(range(10 55)) ylabel(10(5)50, format(%5.0f) labsize(small)) /// 
		xline(13, lpattern(shortdash) lcolor(black) lwidth(medium)) ///
		xlabel(1 `""Jan" "2016" "'  					 ///
				5"May"  9"Sep"   	 ///
				13 `""Jan" "2017" "'					 ///
				17"May" 21"Sep"  ///
				25 `""Jan" "2018" "'					 ///
				29"May"  33"Sep"   ///
				,  format(%tmMon) labsize(small))  ///
		xtitle("", height(4)) ///
		ytitle("Coverage Rate (%)", size(medsmall)) ///
		title("Panel A", size(medsmall) height(0)) ///
		legend(label(1 "Control") label(2 "Treatment") rows(2) pos(0) bplace(neast) order(2 1)) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		saving(1a, replace)

* Panel B

mat drop _all
forval c = 1/36 {
estimates use "$path/Estimates/monthly_cov_effect_4564_051720.ster", number(`c')
mat results = nullmat(results) \ (_b[treatment]*100, 100*(_b[treatment] - 1.96*_se[treatment]) ,100*(_b[treatment] + 1.96*_se[treatment]))
}

clear
set obs 36
gen month = _n
svmat results

rename results1 diff
rename results2 lower 
rename results3 upper

* Figure
twoway (connected diff month, msymbol(o) mcolor(gs1) clcolor(gs1)) ///
	   (rcap upper lower month,  lcolor(gs1) lwidth(vthin)), ///
	    ysc(range(0 2.2)) ylabel(0(0.4)2.0, format(%5.1f) labsize(small)) ///
		xline(13, lpattern(shortdash) lcolor(black) lwidth(medium)) ///
		xlabel(1 `""Jan" "2016" "'  					 ///
				5"May"  9"Sep"   	 ///
				13 `""Jan" "2017" "'					 ///
				17"May" 21"Sep"  ///
				25 `""Jan" "2018" "'					 ///
				29"May"  33"Sep"   ///
				,  format(%tmMon) labsize(small))  ///
		xtitle("", height(4)) ///
		ytitle("Difference in Coverage Rate (p.p.)", size(medsmall)) ///
		title("Panel B", size(medsmall) height(0)) ///
		legend(label(1 "Difference") label(2 "95% CI") ///
					   subtitle("") rows(2) pos(0) bplace(neast)) ///	
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		saving(1b, replace)

		
local note = "{it:Notes}: Panel A displays the shares of middle age adults in the treatment and control groups that were enrolled in any coverage in the specified month. Panel B displays the difference in the share of middle age adults between the treatment and control groups that were enrolled in any coverage in the specified month. Units are percentage points (0-100). Both panels are limited to individuals between the ages of 45 and 64 at the end of 2017 and exclude individuals with full coverage in January through November of 2016. Brackets denote the 95% confidence interval based on standard errors that are clustered by household." 
wordwrap "`note'", l(88)
		
* Combine Panel A and B
* ----------------------
graph combine 1a.gph 1b.gph, col(1) ysize(9) xsize(7) graphregion(margin(l=5 r=5)) ///
	note("`r(text)'",   size(vsmall)  height(17) linegap(0.5)) ///
	graphregion(fcolor(white) color(white) lcolor(white))

graph export "$output/Coverage Effect Over Time (45-64).eps", replace as(eps) 
cap erase 1a.gph 
cap erase 1b.gph 

/* ===============================
Figure A.IV: Coverage Effect by Type of Coverage
==================================*/
clear all
mat drop _all

forval c = 1/36 {
estimates use "$path/Estimates/monthly_medicaid_effect_051720.ster", number(`c')
mat results_med = nullmat(results_med) \ (_b[treatment]*100, 100*(_b[treatment] - 1.96*_se[treatment]) ,100*(_b[treatment] + 1.96*_se[treatment]))
}

forval c = 1/36 {
estimates use "$path/Estimates/monthly_exchange_effect_051720.ster", number(`c')
mat results_exch = nullmat(results_exch) \ (_b[treatment]*100, 100*(_b[treatment] - 1.96*_se[treatment]) ,100*(_b[treatment] + 1.96*_se[treatment]))
}

forval c = 1/36 {
estimates use "$path/Estimates/monthly_esi_effect_051820.ster", number(`c')
mat results_esi = nullmat(results_esi) \ (_b[treatment]*100, 100*(_b[treatment] - 1.96*_se[treatment]) ,100*(_b[treatment] + 1.96*_se[treatment]))
}

set obs 36
gen month = _n

svmat results_med 
svmat results_exch
svmat results_esi

rename results_med1 diff_med
rename results_med2 lower_med 
rename results_med3 upper_med

rename results_exch1 diff_exch
rename results_exch2 lower_exch 
rename results_exch3 upper_exch

rename results_esi1 diff_esi
rename results_esi2 lower_esi 
rename results_esi3 upper_esi

local note = "{it:Notes}: The figure displays the difference in the share of the treatment and control groups enrolled in Medicaid, Exchange and employer-sponsored coverage in the specified month. Units are percentage points (0-100). The figure excludes individuals with full coverage in January through November of 2016. Brackets denote the 95% confidence interval based on standard errors that are clustered by household."
wordwrap "`note'", l(90)	


* Figure
twoway (connected diff_med month , msymbol(o) mcolor(gs1) clcolor(gs1)) ///
	   (rcap upper_med lower_med month,  lcolor(gs1) lwidth(vthin)) ///
	   (connected diff_esi month, msymbol(o) mcolor(gs1) clcolor(gs1) lpattern(shortdash_dot)) ///
	   (rcap upper_esi lower_esi month,  lcolor(gs1) lwidth(vthin)) ///
	   (connected diff_exch month, msymbol(o) mcolor(gs1) clcolor(gs1) lpattern(dash)) ///
	   (rcap upper_exch lower_exch month ,  lcolor(gs1) lwidth(vthin)), ///
		ysc(range(0 1.1)) ylabel(0(0.2)1.0, format(%5.1f) labsize(small)) /// 
		xline(13, lpattern(shortdash) lcolor(black) lwidth(medium)) ///
		xlabel( 1 `""Jan" "2016" "'  					 	///
				3"Mar" 5"May" 7"Jul" 9"Sep" 11"Nov"  	 ///
				13 `""Jan" "2017" "'					 ///
				15"Mar" 17"May" 19"Jul" 21"Sep" 23"Nov"  ///
				25 `""Jan" "2018" "'					 ///
				27"Mar" 29"May" 31"Jul" 33"Sep" 35"Nov"  ///
				,  format(%tmMon) labsize(small))  ///
		xtitle("", height(4)) ///
		ytitle("Difference in Coverage Rate (p.p.)", size(medsmall) height(4)) ///
		legend(label(1 "Medicaid") label(3 "ESI") label(5 "Exchange") label(2 "95% CI") subtitle("") ///
			   order(1 5 3 2) size(small) rows(4) pos(0) bplace(nwest) ) ///
		note("`r(text)'" ,  ///
							size(small) height(15) linegap(0.5)) ///
		graphregion(fcolor(white) color(white) lcolor(white)) 

graph export "$output/Coverage Effect Over Time - Type of Coverage.eps", replace as(eps)

/* ============================================
Figure A.V: Probability of Detecting Mortality Effect by Age Range
=============================================== */


use "$path/Data/power_results_by_age_102619.dta", clear 

* Complier Scalar 2 
twoway (connected power2_01 age, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)), ///
		xtitle("Minimum Age", height(4))		///
		ytitle("Power (%)", height(4) ) ///
		ysize(2) xsize(1.75) ///
		ylabel(5.5(1)7.5) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		title("Percentage Effect: 1%", size(medsmall) ring(1)) ///
		saving(1a, replace)
twoway (connected power2_03 age, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)), ///
		xtitle("Minimum Age", height(4))		///
		ytitle("Power (%)", height(4)) ///
		ysize(2) xsize(1.75) ///
		ylabel(8(4)16) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		title("Percentage Effect: 3%", size(medsmall) ring(1)) ///
		saving(1b, replace)		
twoway (connected power2_05 age, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)), ///
		xtitle("Minimum Age", height(4))		///
		ytitle("Power (%)", height(4) ) ///
		ysize(2) xsize(1.75) ///
		ylabel(15(7.5)30) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		title("Percentage Effect: 5%", size(medsmall) ring(1)) ///
		saving(1c, replace)		
twoway (connected power2_07 age, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)), ///
		xtitle("Minimum Age", height(4))		///
		ytitle("Power (%)", height(4) ) ///
		ysize(2) xsize(1.75) ///
		ylabel(30(10)50) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		title("Percentage Effect: 7%", size(medsmall) ring(1)) ///
		saving(1d, replace)		
		
graph combine 1a.gph 1b.gph 1c.gph 1d.gph, col(2) ysize(10.5) xsize(7.5) graphregion(margin(l=5 r=5)) ///
					graphregion(fcolor(white) color(white) lcolor(white)) ///
					title("Uninsured Complier Scalar: 2", size(medsmall) pos(6) ring(1)) ///
					saving(2, replace)
					
cap erase 1a.gph 
cap erase 1b.gph 
cap erase 1c.gph 
cap erase 1d.gph 



* Complier Scalar 3 
twoway (connected power3_01 age, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)), ///
		xtitle("Minimum Age", height(4))		///
		ytitle("Power (%)", height(4) ) ///
		ysize(2) xsize(1.75) ///
		ylabel(6(0.75)7.5) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		title("Percentage Effect: 1%", size(medsmall) ring(1)) ///
		saving(1a, replace)
twoway (connected power3_03 age, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)), ///
		xtitle("Minimum Age", height(4) )		///
		ytitle("Power (%)", height(4) ) ///
		ysize(2) xsize(1.75) ///
		ylabel(15(5)25) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		title("Percentage Effect: 3%", size(medsmall) ring(1)) ///
		saving(1b, replace)		
twoway (connected power3_05 age, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)), ///
		xtitle("Minimum Age", height(4) )		///
		ytitle("Power (%)", height(4) ) ///
		ysize(2) xsize(1.75) ///
		ylabel(35(10)55) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		title("Percentage Effect: 5%", size(medsmall) ring(1)) ///
		saving(1c, replace)		
twoway (connected power3_07 age, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)), ///
		xtitle("Minimum Age", height(4))		///
		ytitle("Power (%)", height(4)) ///
		ysize(2) xsize(1.75) ///
		ylabel(50(20)90) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		title("Percentage Effect: 7%", size(medsmall) ring(1)) ///
		saving(1d, replace)		
		
graph combine 1a.gph 1b.gph 1c.gph 1d.gph, col(2) ysize(10.5) xsize(7.5) graphregion(margin(l=5 r=5)) ///
					graphregion(fcolor(white) color(white) lcolor(white)) ///
					title("Uninsured Complier Scalar: 3", size(medsmall) pos(6) ring(1)) ///
					saving(3, replace) 
					
cap erase 1a.gph 
cap erase 1b.gph 
cap erase 1c.gph 
cap erase 1d.gph 

*Complier Scalar 4
twoway (connected power4_01 age, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)), ///
		xtitle("Minimum Age", height(4) )		///
		ytitle("Power (%)", height(4)) ///
		ysize(2) xsize(1.75) ///
		ylabel(6(2)10) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		title("Percentage Effect: 1%", size(medsmall) ring(1)) ///
		saving(1a, replace)
twoway (connected power4_03 age, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)), ///
		xtitle("Minimum Age", height(4))		///
		ysize(2) xsize(1.75) ///
		ylabel(25(10)45) ///
		ytitle("Power (%)", height(4) ) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		title("Percentage Effect: 3%", size(medsmall) ring(1)) ///
		saving(1b, replace)		
twoway (connected power4_05 age, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)), ///
		xtitle("Minimum Age", height(4) )		///
		ytitle("Power (%)", height(4) ) ///
		ysize(2) xsize(1.75) ///
		ylabel(50(15)80) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		title("Percentage Effect: 5%", size(medsmall) ring(1)) ///
		saving(1c, replace)		
twoway (connected power4_07 age, msymbol(o) mlcolor(gs1) mfcolor(white) clcolor(gs1) lpattern(dash)), ///
		xtitle("Minimum Age", height(4))		///
		ytitle("Power (%)", height(4)) ///
		ysize(2) xsize(1.75) ///
		ylabel(80(10)100) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		title("Percentage Effect: 7%", size(medsmall) ring(1)) ///
		saving(1d, replace)		
		
graph combine 1a.gph 1b.gph 1c.gph 1d.gph, col(2) ysize(10.5) xsize(7.5) graphregion(margin(l=5 r=5)) ///
					graphregion(fcolor(white) color(white) lcolor(white)) ///
					title("Uninsured Complier Scalar: 4", ring(1) size(medsmall) pos(6)) ///
					saving(4, replace) 
					
					
cap erase 1a.gph 
cap erase 1b.gph 
cap erase 1c.gph 
cap erase 1d.gph 

*** Combine all panels 
graph combine  2.gph 3.gph 4.gph, col(1) ///
			ysize(16) xsize(10) ///
			graphregion(fcolor(white) color(white) lcolor(white)) iscale(*0.73) ///
			note("{it:Notes:} The figure displays the probability of detecting a difference between treatment and control group" ///
				 "mortality at a 5% level of significance, for varying age ranges and mortality effect sizes. Results are " ///
				 "based on simulations with N=1000 random draws of treatment and control populations. Within a figure, " /// 
				 "the x-axis reflects the minimum age included in the analysis; the maximum age is held fixed at 64." ///
				 "The effect of the intervention on coverage for each age range comes from comparing average months of" ///
				 "2017-18 coverage in the treatment and control groups, for individuals with the specified ages that" ///
				 "did not have coverage in each month of 2016. The effect of coverage on mortality is alternatively" /// 
				 "assumed to be a reduction in baseline complier mortality of 1, 3, 5, or 7%. General population mortality" ///
				 "is estimated from population-level mortality rates for 2016 from the Social Security Death Index among" ///
				 "individuals alive at the end of 2015, aggregated into 5-year age bins. The mortality rate for uninsured" ///
				 "compliers, absent insurance, is alternatively assumed to be 2, 3, or 4 times the average mortality" ///
				 "rate for the general population.", linegap(0.5) height(28) size(vsmall)) 
	
graph display, ysize(16) xsize(10)

graph export "$output/Power Analysis - Age Cutoff.eps", replace as(eps)



/*==========================================
Figure A.VI: Effect of Intervention on Additional Middle Age Deaths by Period
==============================================*/

clear all

forval c = 1/6{
estimates use "$path/Estimates/mort_interval_052020.ster", number(`c')
mat estimates = nullmat(estimates) \ (_b[treatment]*100, 100*(_b[treatment] - 1.96*_se[treatment]) ,100*(_b[treatment] + 1.96*_se[treatment]))
}

set obs 6
gen month= 1 + (_n-1)*6
svmat estimates

rename estimates1 percent_dead_diff
rename estimates2 lower
rename estimates3 upper

local note = "{it:Notes}: The figure displays the difference in the mortality rates between middle age adults in the control and treatment groups during the six-month interval extending through the end of the specified month. Units are percentage points (0-100). The analysis is limited to individuals between the ages of 45-64 at the end of 2017 and excludes individuals with full coverage in January through November of 2016. Brackets denote the 95% confidence interval based on standard errors that are clustered by household."
wordwrap "`note'", l(92)	

	
	* Figure
	twoway (scatter percent_dead month , clpattern(1) clcolor(gs1) msymbol(o) mcolor(black)) ///
		   (rcap upper lower month, clpattern(1) lcolor(black) lwidth(thin)), ///
			xline(10, lpattern(shortdash) lcolor(black) lwidth(medium)) ///
			yline(0, lwidth(thin) lcolor(black)) ///
			ysc(range(-0.05 0.05)) ylabel(-0.04(0.02)0.06, format(%5.2f) labsize(small)) ///
			xlabel(1  `" "Jan-Jun" "2016" "'  7  `" "Jul-Dec" "2016" "' ///
				   13 `" "Jan-Jun" "2017" "'  19 `" "Jul-Dec" "2017" "' ///
				   25 `" "Jan-Jun" "2018" "'  31 `" "Jul-Dec" "2018" "' , ///
				   labsize(small)) ///			
			xtitle("", height(0)) ///
			ytitle("Difference in Rate of New Deaths (p.p.)", size(medsmall) height(4)) ///	
			note("`r(text)'") ///
			graphregion(fcolor(white)) ///
	        legend(off)

graph export "$output/Effect of Intervention on Rate of New Deaths.eps", replace as(eps)
	

/*==============
Figure A.VII: Permutation Test for Effect of Intervention on Middle Age Mortality
================*/

use "$path/Data/placebo_dist_110419.dta", clear

rename _pm_3 b_dead1718
rename _pm_4 se_dead1718
replace b_dead1718 = b_dead1718 * 100
replace se_dead1718 = se_dead1718 * 100
gen t_dead1718 = b_dead1718/se_dead1718

*obtain p-values 2017-18
gen extreme1718 = 0
replace extreme1718 = 1 if abs(t_dead1718) >=  2.58
summ extreme1718
local p_val1718 = r(mean)
display `p_val1718'

*Figure
local note = "{it:Notes}: The figure plots the distribution of t-statistics corresponding to the estimated intent-to-treat effect of the intervention on 2017-18 mortality, generated from 1000 random reassignments of the treatment variable across households in the sample population. The vertical line denotes the t-statistic estimated using the actual treatment assignment. The analysis is limited to individuals between the ages of 45-64 at the end of 2017 and excludes individuals with full coverage in January through November of 2016."
wordwrap "`note'", l(90)

local obs =  - 2.58
		histogram t_dead1718, ///
			fraction bin(50) ///
			lcolor(black) fcolor(black)  ///
			xlabel(,labsize(small) format(%12.0f)) ///
			ylabel(,labsize(small) format(%12.2f)) ///
			xtitle("Placebo t-statistics", height(5)) ///
			ytitle("Density", height(4)) ///
			legend(off) ///
			xline(`obs', lcolor(red) ) ///
			text(0.04 -2.6 "Sample" "t-stat", size(small) place(west) just(right)) ///
			note("`r(text)'" , ///
			linegap(0.5) height(20)) ///
			graphregion(fcolor(white))

graph export "$output/Randomization Inference.eps", replace as(eps)



/*==========================================
Figure A.VIII: Effect of Intervention on Middle Age Mortality Among Prior-Year Insured
==============================================*/

use "$path/Data/mortality_by_month_051920.dta", clear

* Keep Placebo Group Only
keep if notall16 == 0

* Calculate rates
forval x = 1/36 {
	replace dead_by_`x' = (dead_by_`x' / dep_tin)*100
}

* Reshape
reshape long dead_by_, i(treatment) j(month) 
rename dead_by_ percent_dead
format dep_tin %12.0f
rename dep_tin count

* Calculate standard error 
gen se = sqrt((percent_dead * (100-percent_dead)) / count)

* Calculate standard error for difference and CI limits
gen upper = percent_dead + 1.96*se
gen lower = percent_dead - 1.96*se

local note = "{it:Notes}: The figure displays the percentage of middle age adults with full prior-year coverage that died during or prior to the specified month. The analysis is limited to individuals between the ages of 45-64 at the end of 2017 that were enrolled in coverage during each of the first 11 months of 2016."
wordwrap "`note'", l(90)	

	* Figure
	twoway (connected percent_dead month if treatment == 1, clpattern(1) clcolor(gs1) msymbol(o) msize(small) mcolor(gs1)) ///
		   (connected percent_dead month if treatment == 0, clpattern(dash) clcolor(black)  msymbol(oh) msize(small) mcolor(gs1)), ///
			xline(13, lpattern(shortdash) lcolor(black) lwidth(medium)) ///
			ysc(range(0 1.5)) ylabel(0(0.2)1.4, format(%5.1f) labsize(small)) ///
			xlabel(1 `" "Jan" "2016" "'  3 `" "Mar" "' ///
				   5 `" "May"  "' 		 7 `" "Jul"  "' ///
				   9 `" "Sep"  "' 		 11 `" "Nov"  "' ///
				   13 `" "Jan" "2017" "' 15 `" "Mar"  "' ///
				   17 `" "May"  "' 		 19 `" "Jul"  "' ///
				   21 `" "Sep" "' 		 23 `" "Nov"  "' ///
				   25 `" "Jan" "2018" "' 27 `" "Mar"  "' ///
				   29 `" "May"  "'		 31 `" "Jul"  "' ///
				   33 `" "Sep" "' 		 35 `" "Nov"  "', ///
					labsize(small)) ///			
			xtitle("", height(0)) ///
			ytitle("Cumulative Mortality (%)", size(medsmall) height(4)) ///
			legend(label(1 "Treatment") label(2 "Control") ///
						   subtitle("") size(small) rows(2) pos(0) bplace(nwest)) ///	
			graphregion(fcolor(white))  ///
			note("`r(text)'")
			
graph export "$output/Cumulative Mortality - Placebo Sample.eps", replace as(eps) 



/* ============================================
Figure A.IX: Coverage and Mortality Effects by Treatment Arm Among Middle Age Adults
=============================================== */

use "$path/Data/arms_collapse_spanish_041620.dta", clear

rename dep_tin count
drop any_covered2017 covered2017 any_covered1718
drop if cov2016_11 == 1
drop if age_45_64 == 0 
egen treatment_arm = group(base_no_sp early_no_sp nonpersonalize_no_sp exemption_info_no_sp base_sp early_sp nonpersonalize_sp exemption_info_sp)

preserve
keep if treatment_arm == 1
rename count count_c
rename covered1718 covered1718_c
rename dead1718 dead1718_c
tempfile control_group
save "`control_group'"
restore 

keep if treatment_arm != 1
merge m:1 cov2016_11 age_45_64 using "`control_group'", nogen

foreach i of var covered1718 dead1718{
gen effect_`i' = `i' - `i'_c 
}

replace effect_dead1718 = effect_dead1718*100

gen effect_dead_se = 100*sqrt(dead1718*(1-dead1718)/count + dead1718_c*(1-dead1718_c)/count_c) ///

gen effect_dead1718_lower = effect_dead1718 - 1.96*effect_dead_se
gen effect_dead1718_upper = effect_dead1718 + 1.96*effect_dead_se

egen sum = sum(count)

gen weight = count/sum 

reg effect_dead1718 effect_covered1718 [pw = weight], robust

local note = "{it:Notes}: The figure plots the estimated effect of the intervention on coverage and mortality, separately for each treatment arm. The coverage effect (x-axis) corresponds to the effect of the intervention on the number of months of coverage enrolled in during 2017-18. The mortality effect (y-axis) corresponds to the effect of the intervention on mortality in 2017-18; units are percentage points (0-100). All analyses are limited to individuals between the ages of 45-64 at the end of 2017 and exclude individuals with full coverage in January through November of 2016. The best linear fit (dashed line), weighted by cell sample size, has a slope of  -0.088, with standard error 0.094." 
wordwrap "`note'", l(93)

twoway 	 (lfit effect_dead1718 effect_covered1718 [pw = weight], lpattern(dash) lcolor(black) lwidth(thin)) ///
		 (scatter effect_dead1718 effect_covered1718 [fw = count], ///
			clpattern(1) clcolor(gs1) msize(vsmall) msymbol(circle) mcolor(black)), ///
			graphregion(fcolor(white) color(white) lcolor(white)) ///
			xtitle("Coverage Effect (Months)", height(4)) ///
			xsc(r(0 0.6)) xlabel(0(0.2)0.6) ///
			ysc(r(-0.14 0.02))  ylabel(-0.12(0.04)0.02) ///
			ytitle("Mortality Effect (p.p.)", height(5)) ///
			legend(off) ///
			note("`r(text)'")
			
graph export "$output/Effect by Treatment Arms.eps", replace as(eps) 


/* ============================================
Figure A.X: Coverage and Mortality Effects by State Among Middle Age Adults
=============================================== */

clear all
mat drop _all
estimates clear

set obs 51
gen state_num = _n

*Coverage Months
forval state = 1/51{
	estimates use "$path/Estimates/cov_state_3_031720", number(`state')
	mat cov = nullmat(cov) \ (_b[treatment] , _se[treatment])
}

svmat cov
rename cov1 months_cov1718_no_controls

*Mortality Effect Results
forval state = 1/51{
	estimates use "$path/Estimates/mort_state_1_031720", number(`state')
	mat mort = nullmat(mort) \ (_b[treatment] , _se[treatment], e(nobs))
}
svmat mort
replace mort1 = mort1*100
replace mort2 = mort2*100

rename mort_11 mort_effect_no_controls
rename mort_13 count

*Weights
egen total_count = sum(count)
gen weights = count/total_count

*Coverage 2017-18 months and Mortality Effect (no controls) 
gen mort_range = (mort_effect_no_controls >= -0.5 & mort_effect_no_controls <= 0.5)

reg mort_effect_no_controls months_cov1718_no_controls [pw = count]
local int = _b[_cons]
local slope = _b[months_cov1718_no_controls]

local note = "{it:Notes}: The figure plots the estimated effect of the intervention on coverage and mortality, separately by state. The coverage effect (x-axis) corresponds to the effect of the intervention on the number of months of coverage enrolled in during 2017-18. The mortality effect (y-axis) corresponds to the effect of the intervention on mortality in 2017-18; units are percentage points (0-100). All analyses are limited to individuals between the ages of 45-64 at the end of 2017 and exclude individuals with full coverage in January through November of 2016.The best linear fit (dashed line) is weighted by population; its slope is -0.055, with standard error 0.120. For ease of presentation, the figure does not plot 4 states with estimated mortality effects greater than 0.5 percentage points or less than -0.5 percentage points, although those states are used in estimating the best linear fit. Collectively, the excluded states represent approximately 1.97 percent of the United States population."
wordwrap "`note'", l(90)

twoway 	   (scatter mort_effect_no_controls months_cov1718_no_controls [fw = count] if mort_range == 1, clpattern(1) clcolor(gs1) msize(tiny) msymbol(o) mcolor(black)) ///
		   (function y = `int' + `slope'*x, range(-0.3 0.6) lpattern(dash) lcolor(black) lwidth(thin)), ///
			graphregion(fcolor(white) color(white) lcolor(white)) ///
	        legend(off) ///
			xtitle("Coverage Effect (Months)", height(4)) ///
			ytitle("Mortality Effect (p.p.)", size(medsmall) height(5)) 	///
			note("`r(text)'")
			
graph export "$output/Effect by State.eps", replace as(eps) 

/* =================================================
Figure A.XI: Cumulative Distribution of Coverage-Months Among Middle Age Adults
====================================================*/

use "$path/Data/count_by_month_4564_1718_061020.dta", clear 
  	
	drop if covered1718 == .
	collapse (sum) count_ind, by(treatment covered1718)	
	
	* Calculate percentage
	qui sum count_ind if treatment == 0, d
	local total = r(sum)
	gen percentage = (count_ind/`total')*100 if treatment == 0
	qui sum count_ind if treatment == 1, d
	local total = r(sum)
	replace percentage = (count_ind/`total')*100 if treatment == 1

	* Calculate cumulative percentage
	bysort treatment: gen cumulative = sum(percentage)
	
*Figure
local note = "{it:Notes}: The figure displays the cumulative distribution function of the number of months of coverage enrolled in during 2017 and 2018 separately for the treatment and control groups. The figure is limited to individuals between the ages of 45-64 at the end of 2017 and excludes individuals with full coverage in January through November of 2016."
wordwrap "`note'", l(90)
	
	* Plot
	twoway (connected cumulative covered1718 if treatment == 1  & covered1718 != 24, ///
			clpattern(l) clcolor(gs1) msymbol(O) mcolor(gs1) msize(small)) ///
			(connected cumulative covered1718 if treatment == 0  & covered1718 != 24, ///
			clpattern(dash) clcolor(gs1) msymbol(oh) mcolor(gs1) msize(small)), ///
			xsc(range(0 23)) xlabel(0(2)23,  labsize(small)) ///
			ysc(range(45 82)) ylabel(50(10)80,labsize(small)) ///
			xtitle("Covered Months", height(4)) ///
			ytitle("Percentage", height(4)) ///
			legend(label(1 "Treatment") label(2 "Control") rows(2) pos(0) bplace(seast)) ///
			title(`cdf`x'', color(black)) ///
			graphregion(fcolor(white) color(white) lcolor(white)) ///
			note("`r(text)'", linegap(0.5) height(15))	
						
graph export "$output/Coverage CDF.eps", replace as(eps)

/*==============
Figure A.XII: Cumulative Distribution of Coverage-Months Among Middle Age Adults by Demographic Group
================*/

* Panel A and B: by gender 
use "$path/Data/cov1718_counts_by_gender.dta", clear //gender

*Clean
drop if covered1718 == .
collapse (sum) count_ind, by(female treatment covered1718)
sort female treatment covered1718

* Create percentages
gen N = .
gen pct_covered = .
forvalues x = 0/1 {
	forvalues f = 0/1 {
		qui sum count_ind if treatment == `x' & female == `f'
		replace N = r(sum) if treatment == `x' & female == `f'
		replace pct_covered = count_ind / N * 100 if treatment == `x' & female == `f'
	}
}
					
* Calculate cumulative percentage
bysort female treatment: gen cumulative = sum(pct_covered)
	
* Plot
twoway (connected cumulative covered1718 if treatment == 1  & covered1718 != 24 & female == 0, ///
		clpattern(l) clcolor(gs1) msymbol(O) mcolor(gs1) msize(small)) ///
	   (connected cumulative covered1718 if treatment == 0  & covered1718 != 24 & female == 0, ///
		clpattern(dash) clcolor(gs1) msymbol(oh) mcolor(gs1) msize(small)), ///
		xsc(range(0 23)) xlabel(0(2)23,  labsize(small)) ///
		ysc(range(45 95)) ylabel(50(10)90,labsize(small)) ///
		xtitle("Covered Months", height(3)) ///
		ytitle("Percentage", height(3)) ///
		legend(label(1 "Treatment") label(2 "Control") size(small) rows(1) pos(0) bplace(seast)) ///
		title(Men, color(black) size(medium)) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
	saving(a, replace)


* Plot
twoway (connected cumulative covered1718 if treatment == 1  & covered1718 != 24 & female == 1, ///
		clpattern(l) clcolor(gs1) msymbol(O) mcolor(gs1) msize(small)) ///
		(connected cumulative covered1718 if treatment == 0 & covered1718 != 24 & female == 1, ///
		clpattern(dash) clcolor(gs1) msymbol(oh) mcolor(gs1) msize(small)), ///
		xsc(range(0 23)) xlabel(0(2)23,  labsize(small)) ///
		ysc(range(45 95)) ylabel(50(10)90,labsize(small)) ///
		xtitle("Covered Months", height(3)) ///
		ytitle("Percentage", height(3)) ///
		title(Women, color(black) size(medium)) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		saving(b, replace)

*** Panel C and D: by marriage
use "$path/Data/cov1718_counts_by_married.dta", clear //married

*Clean
drop if covered1718 == .
collapse (sum) count_ind, by(married treatment covered1718)
sort married treatment covered1718

* Create percentages
gen N = .
gen pct_covered = .
forvalues x = 0/1 {
	forvalues m = 0/1 {
		qui sum count_ind if treatment == `x' & married == `m'
		replace N = r(sum) if treatment == `x' & married == `m'
		replace pct_covered = count_ind / N * 100 if treatment == `x' & married == `m'
	}
}
					
* Calculate cumulative percentage
bysort married treatment: gen cumulative = sum(pct_covered)
	
* Plot
twoway (connected cumulative covered1718 if treatment == 1  & covered1718 != 24 & married == 0, ///
		clpattern(l) clcolor(gs1) msymbol(O) mcolor(gs1) msize(small)) ///
	   (connected cumulative covered1718 if treatment == 0  & covered1718 != 24 & married == 0, ///
		clpattern(dash) clcolor(gs1) msymbol(oh) mcolor(gs1) msize(small)), ///
		xsc(range(0 23)) xlabel(0(2)23,  labsize(small)) ///
		ysc(range(45 95)) ylabel(50(10)90,labsize(small)) ///
		xtitle("Covered Months", height(3)) ///
		ytitle("Percentage", height(3)) ///
		title(Married, color(black) size(medium)) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		saving(c, replace)


* Plot
twoway (connected cumulative covered1718 if treatment == 1  & covered1718 != 24 & married == 1, ///
		clpattern(l) clcolor(gs1) msymbol(O) mcolor(gs1) msize(small)) ///
		(connected cumulative covered1718 if treatment == 0 & covered1718 != 24 & married == 1, ///
		clpattern(dash) clcolor(gs1) msymbol(oh) mcolor(gs1) msize(small)), ///
		xsc(range(0 23)) xlabel(0(2)23,  labsize(small)) ///
		ysc(range(45 95)) ylabel(50(10)90,labsize(small)) ///
		xtitle("Covered Months", height(3)) ///
		ytitle("Percentage", height(3)) ///
		title(Not Married, color(black) size(medium)) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		saving(d, replace)



*** Panel E and F: by income
use "$path/Data/cov1718_counts_by_income.dta", clear //income

*Clean
drop if covered1718 == .
collapse (sum) count_ind, by(above_med_inc treatment covered1718)
sort above_med_inc treatment covered1718

* Create percentages
gen N = .
gen pct_covered = .
forvalues x = 0/1 {
	forvalues a = 0/1 {
		qui sum count_ind if treatment == `x' & above_med_inc == `a'
		replace N = r(sum) if treatment == `x' & above_med_inc == `a'
		replace pct_covered = count_ind / N * 100 if treatment == `x' & above_med_inc == `a'
	}
}
					
* Calculate cumulative percentage
bysort above_med_inc treatment: gen cumulative = sum(pct_covered)

* Plot
twoway (connected cumulative covered1718 if treatment == 1  & covered1718 != 24 & above_med_inc == 1, ///
		clpattern(l) clcolor(gs1) msymbol(O) mcolor(gs1) msize(small)) ///
		(connected cumulative covered1718 if treatment == 0 & covered1718 != 24 & above_med_inc == 1, ///
		clpattern(dash) clcolor(gs1) msymbol(oh) mcolor(gs1) msize(small)), ///
		xsc(range(0 23)) xlabel(0(2)23,  labsize(small)) ///
		ysc(range(45 95)) ylabel(50(10)90,labsize(small)) ///
		xtitle("Covered Months", height(3)) ///
		ytitle("Percentage", height(3)) ///
		title(High Income, color(black) size(medium)) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		saving(e, replace)
	
* Plot
twoway (connected cumulative covered1718 if treatment == 1  & covered1718 != 24 & above_med_inc == 0, ///
		clpattern(l) clcolor(gs1) msymbol(O) mcolor(gs1) msize(small)) ///
	   (connected cumulative covered1718 if treatment == 0  & covered1718 != 24 & above_med_inc == 0, ///
		clpattern(dash) clcolor(gs1) msymbol(oh) mcolor(gs1) msize(small)), ///
		xsc(range(0 23)) xlabel(0(2)23,  labsize(small)) ///
		ysc(range(45 95)) ylabel(50(10)90,labsize(small)) ///
		xtitle("Covered Months", height(3)) ///
		ytitle("Percentage", height(3)) ///
		title(Low Income, color(black) size(medium)) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		saving(f, replace)


*Figure
local note = "{it:Notes}: Each figure displays the cumulative distribution function of the number of months of coverage enrolled in during 2017 and 2018 among individuals in the specified demographic group. Individuals are categorized as high or low income based on whether their 2015 adjusted gross income was above or below the sample median. The figure is limited to individuals between the ages of 45-64 at the end of 2017 and excludes individuals with full coverage in January through November of 2016."
wordwrap "`note'", l(93)

*** Combine all six panels 
grc1leg  a.gph b.gph c.gph ///
			  d.gph e.gph f.gph, col(2) ///
			 ysize(10.4) xsize(8) ///
	note("`r(text)'", ///
					size(vsmall)  height(15) linegap(0.5)) ///
	legendfrom(a.gph) ///
	position(6) ///
	graphregion(fcolor(white) color(white) lcolor(white)) 

graph display, ysize(15) xsize(12)
				
graph export "$output/Coverage CDF by Subgroup.eps", replace as(eps) 

cap erase a.gph
cap erase b.gph
cap erase c.gph
cap erase d.gph
cap erase e.gph
cap erase f.gph


/* ============================================
Figure A.XIII: Distribution of Coverage-Months Induced by the Intervention Among Middle Age Adults
=============================================== */

clear all
mat drop _all
estimates clear

forval c = 1/24 {
	estimates use "$path/Estimates/extra_months_4564_110519.ster", number(`c')
	eststo c`c'
	mat weights = nullmat(weights) \ (_b[treatment])
}
svmat weights

qui sum weights1 // turn into percentages
replace weights1 = (weights1 / r(sum))*100

gen month = _n

*Figure
local note = "{it:Notes}: The figure displays the share of coverage-months induced by the intervention that represent the specified month of coverage for an individual during 2017-18. For example, the bar for month 3 indicates that approximately 5% of the additional months of coverage induced by the intervention represent the third month of coverage that an individual was enrolled in during 2017-18. The bar for a particular month is estimated by comparing the share of the treatment group that enrolled in at least the specified number of months of coverage during 2017-18 with the corresponding share of the control group that did so. Repeating this calculation for each month yields each new month's share of all coverage-months induced by the intervention. The figure is limited to individuals between the ages of 45-64 at the end of 2017 and excludes individuals with full coverage in January through November of 2016."
wordwrap "`note'", l(90)

* Plot
graph bar weights1, over(month) intensity(inten70) bar(1,lcolor(black)) ///
		ysc(range(0 6.5)) ylabel(0(2)6.5,labsize(small) format(%12.0f)) ///
		b1title("New Month Number", height(3) size(medsmall)) ///
		ytitle("Share of All New Months (%)", size(medsmall) height(4)) ///
		legend(off) ///
		note("`r(text)'", ///
			 linegap(0.5) height(35)) ///
			 graphregion(fcolor(white) color(white) lcolor(white)) 

graph export "$output/Weight Estimates by Month.eps", replace as(eps)


/* ============================================
Figure A.XIV: Distribution of Coverage-Months Induced by Treatment Arm Among Middle Age Adults
=============================================== */

clear all
mat drop _all
estimates clear

set obs 24
gen month = _n

local arm early nonpersonalize exemption_info base 
			
foreach i in `arm'{
forval c = 1/24 {
	estimates use "$path/Estimates/`i'_months_060220.ster", number(`c')
	eststo c`c'
	mat `i'_weights = nullmat(`i'_weights) \ (_b[treatment])
}
svmat `i'_weights
rename `i'_weights1 `i'_weights
qui sum `i'_weights 
replace `i'_weights = 100*(`i'_weights/ r(sum))
}



* Plot
graph bar base_weights, over(month, label(labsize(tiny)))  intensity(inten70) bar(1,lcolor(black)) ///
		ysc(range(0 6.5)) ylabel(0(2)6.5,labsize(small) format(%12.0f)) ///
		b1title("New Month Number", height(3) size(small)) ///
		ytitle("Share of All New Months (%)", size(small) height(4)) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
		saving(a, replace) ///
		title(Base, size(medsmall))
			 
* Plot
graph bar early_weights, over(month, label(labsize(tiny)))  intensity(inten70) bar(1,lcolor(black)) ///
		ysc(range(0 6.5)) ylabel(0(2)6.5,labsize(small) format(%12.0f)) ///
		b1title("New Month Number", height(3) size(small)) ///
		ytitle("Share of All New Months (%)", size(small) height(4)) ///
		legend(off) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
			 saving(b, replace) ///
		title(Early, size(medsmall))
			 
* Plot
graph bar nonpersonalize_weights, over(month, label(labsize(tiny)))  intensity(inten70) bar(1,lcolor(black)) ///
		ysc(range(0 6.5)) ylabel(0(2)6.5,labsize(small) format(%12.0f)) ///
		b1title("New Month Number", height(3) size(small)) ///
		ytitle("Share of All New Months (%)", size(small) height(4)) ///
		legend(off) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
			 saving(c, replace)  ///
		title(Non-Personalized, size(medsmall))
			 
* Plot
graph bar exemption_info_weights, over(month, label(labsize(tiny))) intensity(inten70) bar(1,lcolor(black)) ///
		ysc(range(0 6.5)) ylabel(0(2)6.5,labsize(small) format(%12.0f)) ///
		b1title("New Month Number", height(3) size(small)) ///
		ytitle("Share of All New Months (%)", size(small) height(4)) ///
		legend(off) ///
		graphregion(fcolor(white) color(white) lcolor(white)) ///
			 saving(d, replace)  ///
		title(Exemption Info, size(medsmall)) 

*Figure
local note = "{it:Notes}: The figure displays the share of coverage-months induced by the specified treatment arm that represent the specified month of coverage for an individual during 2017-18. The bar for a particular month is estimated by comparing the share of the specified treatment arm that enrolled in at least the specified number of months of coverage during 2017-18 with the corresponding share of the control group that did so. Repeating this calculation for each month yields each new month's share of all coverage-months induced by the intervention. The figure is limited to individuals between the ages of 45-64 at the end of 2017 and excludes individuals with full coverage in January through November of 2016.  "
wordwrap "`note'", l(100)			 

*** Combine all panels 
graph combine  a.gph b.gph c.gph d.gph, col(2) ///
			graphregion(fcolor(white) color(white) lcolor(white))  ///
			note("`r(text)'",  size(vsmall))  ///
			ysize(8) xsize(9) graphregion(margin(l=5 r=5))
	

graph export "$output/Weight Estimates by Month - Four Arms.eps", replace as(eps)

cap erase a.gph
cap erase b.gph
cap erase c.gph
cap erase d.gph


